Numerical study of the random field Ising model at zero and positive temperature 
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In this paper the three dimensional random field Ising model is studied at both zero temperature 
and positive temperature. Critical exponents are extracted at zero temperature by finite size scaling 
analysis of large discontinuities in the bond energy. The heat capacity exponent a is found to be 
near zero. The ground states are determined for a range of external field and disorder strength 
near the zero temperature critical point and the scaling of ground state tilings of the field-disorder 
plane is discussed. At positive temperature the specific heat and the susceptibility are obtained 
using the Wang-Landau algorithm. It is found that sharp peaks are present in these physical 
quantities for some realizations of systems sized 16^ and larger. These sharp peaks result from 
flipping large domains and correspond to large discontinuities in ground state bond energies. Finally, 
zero temperature and positive temperature spin configurations near the critical line are found to be 
highly correlated suggesting a strong version of the zero temperature fixed point hypothesis. 

PACS numbers: 75.10.Nr, 05.70.Fh, 75.10.Hk 



I. INTRODUCTION 

The random field Ising model (RFIM) is among the 
simplest non-trivial spin models with quenched disorder. 
It has been intensively studied theoretically, experimen- 
tally, and in computer simulations during the last thirty 
years but is still not well understood. Following the 
seminal discussion of Imry and Mai it has been proved 
that the RFIM has an ordered phase at low temperature 
and weak disorder when the dimension is greater than 
^.^^2.3.4^ It is generally believed that the transition from 
the ordered phase to the disordered phase of the RFIM is 
continuous and is controlled by a zero temperature fixed 
point^'*''^. Since random field fluctuations dominate over 
thermal fluctuations at the transition, the hyperscaling 
relation is modified as {d — 6)1^ = 2 — a, where 9 is the 
violation of hyperscaling exponent^. 

The phase diagram of the RFIM is sketched in Fig. ^ 
Phase transitions can occur from the ferromagnetic phase 
(F) to the paramagnetic phase (P) at either zero temper- 
ature as a function of disorder strength A at A = Ac, 
or as a function of temperature T if disorder is fixed at 
Aq < Ac. According to the zero temperature fixed point 
hypothesis, the zero temperature transition and the pos- 
itive temperature transitions belong to the same univer- 
sality class. In this paper we use numerical methods to 
study both kinds of transitions and connections between 
them. One of our primary results is that, for each real- 
ization of disorder, there is a strong correlation between 
ground state configurations near Ac and thermal states 
near Tc for Aq < Ac. 

Currently, there is no controlled rcnormalization group 
analysis of the RFIM phase transition and Monte Carlo 
simulations of the RFIM at positive temperaturc^iiSiiiiiS 
are limited to small systems because of very long equili- 
bration times^. According to the zero temperature fixed 




FIG. 1: Phase diagram of the RFIM. We study two types of 
phase transitions going from the ferromagnetic phase (F) to 
the paramagnetic phase (P): The zero temperature transition 
(open arrow) occurs at T = and A = Ac, and positive 
temperature transitions (solid arrow) occurs at a fix disorder 
A = Ao < Ac and T > 0. 



point hypothesis, many properties of the RFIM phase 
transition, including the values of critical exponents, can 
be determined by studying the RFIM at zero tempera- 
ture. The ground state of the RFIM can be found in 
polynomial time^'' by efficient combinatorial algorithms 
so that zero temperature simulations are much faster and 
allow for much larger system sizes than positive tempera- 
ture simulations. Critical exponents have been obtained 
from zero temperature studie9i^'iS*i& that are mostly con- 
sistent with the scaling theories^'S'i, series methodsii and 
real space rcnormalization group approachesiS^iSiSS. 
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Much work has been done in determining the critical 
exponents, and the values of many exponents are well 
established. However, the value of the heat capacity ex- 
ponent a is still controversial. A recent zero tempera- 
ture study by Hartmann and Young^** found a « —0.6 
for the three dimensional Gaussian RFIM. The modi- 
fied hyperscaling relation, however, predicts that a = 
2 — {d — 0)v K, 0, given the well-accepted values 6 k, 1.5 
and V K, 1.1 — 1.4. Therefore the quite negative value 
found in Ref. ^3 is inconsistent with the modified hyper- 
scaling relation. Some older work at zero temperaturaSi 
and Monte Carlo simulations^^ also found a quite neg- 
ative. On the other hand, Middleton and Fisher-^ also 
studied the three dimensional Gaussian RFIM at T = 
and found a « —0.1, in agreement with the modified 
hyperscaling relation. 

Dukovski and Machtai^ studied the ground states of 
the RFIM in the presence of an external field H . They 
located a "finite size critical point" for each realization 
of disorder, identified as the point of degeneracy of three 
ground states in the the H — A plane with the largest 
discontinuity in magnetization. They extracted critical 
exponents via finite size scaling of the discontinuities at 
that point. The reason to focus on the finite size crit- 
ical point was that this point can be regarded as the 
most singular point on the H — A plane, and working at 
this point may reduce the influence of the regular part 
of physical quantities. The value of heat capacity expo- 
nent they found was a « 0, however, their results were 
less accurate than those of Refs. IT^ and ITsI because of the 
large amount of computational work needed to locate the 
finite size critical point. 

The work reported in this paper combines both zero 
temperature and positive temperature studies. The zero 
temperature studies extend the work of Ref. in two di- 
rections. First, for each realization of disorder we study 
points along the H — line with large discontinuities 
in bond energy or magnetization to determine the criti- 
cal exponents. Finding discontinuities along the H = 
line requires much less computational work than find- 
ing the finite size critical point while still adhering to the 
idea introduced in Ref. Ha of extracting critical exponents 
from the large discontinuities in each realization of disor- 
der. We also find ground state spin configurations near 
these large discontinuities and compare them to thermal 
states near positive temperature critical points. Second, 
we study the full set of ground states of the RFIM in the 
critical region of the H ~ A plane and discuss the prop- 
erties of the resulting ground state tilings of this plane. 

Since conventional Monte Carlo methods are not ef- 
ficient for the study of the RFIM, we apply the Wang- 
Landau algorithm^^ to the RFIM, which enables us to ob- 
tain the specific heat and the susceptibility over a broad 
range of temperature with system size up to 32"^. We find 
that some realizations display sharp peaks in the specific 
heat and susceptibility. Inspired by the zero temperature 
fixed point hypothesis, we relate these sharp peaks to 
the large discontinuities at zero temperature. We further 



study the thermal states (average spin configurations) 
near the transition using the Metropolis algorithm and 
compare them to the ground states near the zero temper- 
ature transition. Some of this work has been previously 
announced in Ref. |2^ 

In this paper we consider the three dimensional RFIM 
with Gaussian random fields described by the Hamilto- 
nian, 

H = SiSj - A^^hiSi ~ H^^Si, (1) 

where H is the uniform external field, {i,j) indicates a 
sum over all nearest neighbor sites i and j on a simple 
cubic lattice of linear size L with periodic boundary con- 
ditions. The random fields hi are Gaussian random vari- 
ables with mean zero and standard deviation one and the 
strength of disorder is A. The normalized fields {hi} de- 
fine a realization of disorder and, for a given realization 
of disorder we explore spin configurations and physical 
properties as a function of H, T and A. Some of the 
physical quantities of interest include the magnetization 
m, defined as 

(2) 

i 

and the bond energy e, 

e = --^^s,Sj-. (3) 

In the next section we discuss the scaling properties 
of large discontinuities in the bond energy at zero tem- 
perature and use numerical results for these discontinu- 
ities to extract critical exponents and the critical disorder 
strength. In Sec. lIIII we obtain ground state portraits for 
the RFIM and discuss their scaling properties. Section 
IIVI presents results of positive temperature simulations 
and, in Sec. we discuss correlations between ground 
states and thermal states. The paper closes with a sum- 
mary and discussion. 

II. CRITICAL EXPONENTS AT ZERO 
TEMPERATURE 

At zero temperature, the problem of finding the ground 
state of the RFIM can be mapped to the MAX-FLOW 
problem in graph theory, which is solvable in polynomial 
timei^. We use a modified version of the push-relabel 
algorithm^ to calculate the ground state o^'^i^^ . 

The specific heat at T = is definedi^ as 



where e is the bond energy defined in Eq. ((2Jl and the 
square brackets denote averaging over disorder realiza- 
tions. At zero temperature, for each realization of nor- 
malized random fields the bond energy changes discon- 
tinuously as a function of the strength of disorder A. An 
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example of a single realization is shown in Fig. |^ and 
illustrates the point that the sizes of the discontinuities 
vary widely. In Refs. ^3 and all bond energy jumps 
are included in the calculation of the specific heat expo- 
nent. However, small jumps arc presumably part of the 
analytic background rather than the singular behavior. 
The smallest bond energy jump, for example, is 4 for all 
system sizes. We therefore analyze large jumps in the 
bond energy to focus on the critical singularity. From a 
sample of N{L) disorder realizations of systems of size 
L, let [Sei] be the average over the N{L) largest bond 
energy jumps and let [Ai] be the average of the disorder 
strength at these jumps. Each realization of disorder typ- 
ically contributes one bond energy jump and one disorder 
strength to these averages though some realizations con- 
tribute nothing and some several values. Figure |3| shows 
the specific heat decomposed into two components, com- 
ponent (a) is due to the largest N{L) jumps while com- 
ponent (b) arises from all other jumps. One can see that 
the large jumps make a significant contribution to the full 
specific heat and we will use this component to extract 
critical exponents. The finite size scaling of the specific 
heat is expected to obey 

C ^ L°'/''C{{A- Ac)L^^'') (5) 

Though the two components of the specific heat shown 
in Fig. 121 obviously behave differently from one another, 
our primary assumption is that the finite-size scaling of 
the full specific heat also applies to the component of the 
specific heat from the largest bond energy jumps. Indeed 
we believe that the large jumps provide better data to ob- 
tain critical exponents from small systems than the full 
specific heat because this component is undiluted by the 
analytic background. We discuss more detailed scaling 
assumptions about large jumps later in this section. Note 
that the peak height for component (a) barely changes 
with system size suggesting that the specific heat expo- 
nent a is near zero. 

Integrating Eq. |SJ) as applied to the component from 
the largest discontinuities, we obtain a finite size scaling 
ansatz for the large jumps, 

[,5ei] - (6) 

where a is the specific heat exponent and v is the corre- 
lation length exponent. Table gives the average size 
of the large bond energy jumps as a function of sys- 
tem size L. A fit of the form given in Eq. (0 yields 
(1 — oi)lv = 0.842 ± 0.003 with goodness of fit parameter 
Q K. 0.7 (Q = r(d/2, xV2) with d the number of degrees 
of freedom and F the incomplete gamma function) . 

The displacement of the average position of the large 
jumps from the infinite volume limit and the standard 
deviation of the positions of the large jumps are each 
measures of the width of the critical region and, following 
Ref. we assume that they satisfy the finite size scaling 
relations, 

[Ai] - A, ^ L-^l\ (7) 
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FIG. 2: Bond energy as a function of disorder strength A for 
a single realization of disorder (seed 1003). Numbers 1 and 2 
indicate the two biggest jumps in the bond energy. 
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FIG. 3: Two components of the specific heat, (a) The contri- 
bution to the specific heat from the largest N[L) bond energy 
jumps, where N{L) is the number of disorder realizations for 
each system size L. (b) The contribution from all other bond 
energy jumps. 

and 

V[(Ai - [Ai])2] ^ L-'/r (8) 

where v is the correlation length exponent and Ac is the 
infinite size critical disorder strength. Table Ogives the 
standard deviation of the position of the largest jump 
and a fit to Eq. 10 yields l/i/ = 0.79 ± 0.01 with Q w 
0.4. Finally, Table Ogives [Ai] and, using the previously 
obtained value, l/i/ = 0.79 a fit to Eq. yields Ac = 
2.280 ±0.003 with Q « 0.2. 

The second and third largest jumps also presumably 
refiect the critical singularity. We repeated the foregoing 
calculations for the largest kN{L) jumps, where our data 
allow us to go up to fc = 3. The results are listed in Table 
im We arrive at the following best estimates of the critical 
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TABLE I: Data from ground state simulations for the average 
largest jump, and the average and standard deviation of the 
disorder strength at the largest jumps as a function of system 
size L. 



L 




[(Ai - [Ax])^]^/^ 


[Ai] 


32 


0.1208(4) 


0.0904(5) 


2.4915(10) 


48 


0.0857(5) 


0.0666(6) 


2.4326(10) 


64 


0.0675(5) 


0.0524(5) 


2.4017(10) 


96 


0.0480(6) 


0.0385(7) 


2.3709(10) 


TABLE 


II: Critical exponents extracted from the largest 


kN{L) jumps, where N{L) is the number of realizations for 


system size L. Errors are 


purely statistical. 




k 


{l-a)/u 


1/^ 


A, 


1 


0.841(4) 


0.79(1) 


2.282(2) 


2 


0.842(4) 


0.80(1) 


2.282(2) 


3 


0.844(3) 


0.81(1) 


2.283(1) 



exponent and the infinite volume critical disorder, 

i—^ = 0.842 ± 0.004, J/ =1.25 ±0.02 

Ac 2.282 ±0.002, a = -0.05 ± 0.02. (9) 

where the error bars include statistical errors from all 
three k values. 

Our values of the (l — a)/!/ and Ac are consistent with 
some previous calculations and a is found to be near 
zero, which is in agreement with Ref. ITsL But the value 
of v we have calculated is smaller than recent results 
quoted in Refs. ^3 and In Table IIIII our calculated 
values of the exponents are listed in comparison with 
some recent work. Our values of (1 — a)/^ and 1/v gives 
(2 — a)lv K. 1.64. Applying the modified hyperscaling 
relation and the inequality Q > d/2 — (3/v one has 
j3/v > 0.14, which is inconsistent with other work. We 
believe that our value of (1 — a)/v is more reliable than 
our value of 1/v. The fit for 1/v starts from size L = 32 
and would be quite poor if the L = 16 data were included 
suggesting significant finite size correction. On the other 
hand, if the L = 16 data were included, the fit would still 
be good for (1 — a)lv and there would be no change in 
the resulting value. 

Next we take a closer look at the finite size scaling 
properties of the distribution of discontinuities in the 
bond energy. Bond energy jumps result from flipping 
domains. Considering N{L) {N ^ 1) realizations with 
system size L, we assume that there is a renormaliza- 
tion group transformation mapping them to N{L') real- 
izations with system size L', and the flipped domains are 
transformed such that the average bond energy jump [^e] 
and the average disorder where jumps occur [A] conform 
to the already known scaling relations, Eqs. © and (Tjl. 
In order to exclude small jump that do not scale prop- 
erly, we introduce a lower cut-off (5emin for bond energy 



TABLE III: A summary of recent estimates of Ac, ly, {l — a)/v 
and Q, either calculated by ground state (GS) or Monte Carlo 
(MC) simulations. 



Ref. 


Ac V 




a. 


method 


This work 2.282(2) 1.25(2) 


0.842(4) 


-0.05(2) 


GS 




2.29(2) 1.1(1) 


0.80(3) 


0.12 


GS 


ii 


2.28(1) 1.36(1) 


1.20 


-0.63(7) 


GS 




2.270(4) 1.37(9) 


0.82(2)" 


-0.12(12) 


GS 




2.26(1) 1.22(6) 






GS 


2^ 


2.29(4) 1.19(8) 






MC 


21 


2.37(5) 1.0(1) 


1.55 


-0.55(20) 


MC 



"This value was calculated from scaling of the bond energy. They 
found (1 — a)/v = 0.74(2) by relating it to the fractal dimension of 
the surface of spin clusters. 



jumps, which should scale the same way as [5e], 

<5e„,in ~ id-")/". (10) 

The total number of bond energy jumps larger than 
the scaled lower cut-off is independent of the system 
size, since the jumps in systems with different sizes 
are connected by the renormalization group transforma- 
tion. Defining scaled variables u — SeL^^^^"'^^" and 
V ~ {A — Ac)L~^^'^ , the number of jumps occurring in 
a small neighborhood of {u, v) should also be invariant 
for different system sizes. It then follows that the prob- 
ability P{5e, A) of having a bond energy jump with size 
(5e, 5e > Se^i^ and position A is proportional to a given 
normalized probability distribution function P{u,v), 

P{Se, A) cx P((5ei(i-")/^ (A - Ae)^^/'^). (11) 

The normalization of P{Se, A) then gives 

P{6e, A) = L(2~"'/'^P((5ei(i""'/^ (A - Ac)L^^''). (12) 

Letting a = ^emin^^^""'^'' and integrating gives the scal- 
ing of the specific heat due to big jumps Cf,, 

Cb = / 6eP{6e,A)d6e 

"'<5e,„i„ 

^ L"/'^C((A- Ac)Ll/^a), (13) 

where C(v,a) — J^uP{u,v) du is some scaling function. 
By integrating Eq. H12|l over Se we derive the probabil- 
ity distribution of the disorder strength where big bond 
energy jumps occur, 

Pfc(A)-Ll/^Q((A-Ac)Ll/^a), (14) 

where Q{v, a) — P{u, v) du. We then recover the scal- 
ing of the average disorder strength given in Eq. lO The 
setup of a scaled lower cut-off should be equivalent to 
picking the kN{L) largest jumps for each size L, where k 
is some fixed number, because the total number of jumps 
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TABLE IV: Number of bond energy jumps larger than the 
scaled lower cut-off Semm- The cut-off satisfies that a — 
SeininL^^~°'^^'^ is a constant. 







L = 16 


L = 32 


L = 48 


L = 64 


L = 96 


a = 


0.05 


69.3(1) 


75.1(1) 








a = 


0.1 


29.03(6) 


30.89(6) 








a = 


0.6 


3.16(1) 


3.04(1) 


3.01(1) 


2.99(2) 


3.00(3) 


a — 


1 


1.585(5) 


1.533(7) 


1.53(1) 


1.52(1) 


1.54(2) 




FIG. 4: Data collapse of the specific heat. The cut-off is 
a = ^emini'^""'^" = 1.0. System sizes range from 16'' to 96^. 
The inset shows data collapse of the specific heat for system 
size 16^ and 32^ and the cut-off is (Jemini'^""''''' = 0.1. 



is invariant under the renormalization group transfor- 
mation. We test this hypothesis by fixing the constant 
a = SeminL^^~°'^^'^ and count how many jumps larger 
than the lower cut-off there are for each system size. The 
result is hsted in Table Hvl One can see that if the scaled 
lower cut-off is not too small, the number of bond energy 
jumps larger than the cut-off goes to a constant indepen- 
dent of the system size L. 

Figure^illustrates the data collapse of the specific heat 
predicted by Eq. I^Sj) for system sizes 32^, 48^, 64^ and 
96'^. In conventional data collapse the scaling function 
C{x) behaves like x^" as x ^ oo, but in Fig. 01 the tail 
of the curve decays faster than a power law. The main 
plot in Fig.2]has the bond energy jump cut-off set as a = 
'5e„iini(^""^/'' = 1. Tabic HVI shows that for a = 1 only 
a few largest jumps per realization contribute to Cb and, 
since these jumps are concentrated near the critical point, 
we do not expect Cf, ~ (A — Ac)^". However, the tail 
should approach x"", if the number of jumps included is 
increased, or equivalently, a is reduced. The inset in Fig. 
^ shows Cb with a smaller cut-off, a = SeminL^^""^^" — 
0.1. With this cut-off data is available for 16'^ and 32'^ 
systems only. The inset illustrates that, as the cut-off 
is lowered, the tail of the scaling function expands and 
presumably approaches the asymptotic shape. 



III. GROUND STATES PICTURES AND 
SCALING RELATIONS 

The tiling of the H — A plane by ground states is the 
subject of this section. To study this tiling, we find all 
ground states within a certain range of disorder A and 
external field H in the critical region. Since the Hamil- 
tonian of the RFIM is linear with respect to the external 
field H and the strength of disorder A, each spin config- 
uration is represented by a plane in the H — A — Ti. co- 
ordinate system. Ground states are spin configurations 
that are locally lowest and the set of all ground states 
form a convex surface in this coordinate system. Spin 
configurations are ground states within regions of H and 
A bounded by neighboring ground state planes so that 
a given spin configuration is the ground state within a 
polygonal region. At boundaries of these polygons, and 
intersection points of boundaries, ground states are de- 
generate. We are particularly interested in the degener- 
ate points that are common points of three ground states. 
We call these "triple points." 

The structure of the ground state energy surface can 
be visualized by projecting it onto the H — A plane where 
it becomes a tiling of the plane by polygons. The com- 
putational method for finding this tiling is closely related 
to the method developed in Ref. In that work, the 
first order line, which is a set of boundaries that sepa- 
rates the two ordered phases with positive and negative 
magnetization, respectively, was followed and the "finite 
size critical point" was identified by finding the triple 
point having the largest discontinuity in magnetization. 
The finite size critical point was regarded as the most 
singular point, and critical exponents were extracted via 
finite size scaling of magnetization and bond energy dis- 
continuities at the point. In this paper we use a method 
similar to the techniques used in Ref. |2^ and to map 
out all the ground states for any given realizations within 
a certain region on the H — A plane near the finite size 
critical point. 

Our algorithm performs a breadth-first search of 
ground states. The starting point of the search is the 
finite size critical point located by the algorithm of Ref. 

For each point where more than two ground states 
are degenerate we already have the ground states around 
the point and the coexistence lines separating them (one 
locates the point by finding ground states around it). We 
then follow the lines and search for the next adjacent de- 
generate point using the following method. Starting from 
the given degenerate point p, we extend the coexistence 
line separating ground states Pi and P2 with some pre- 
selected step size until we meet a point go on which the 
ground state Qq is different from both Pi and P2- The 
actual adjacent degenerate point is typically passed over, 
because of the the fixed step size is too large. We then 
locate the intersection point of Pi , P2 and Qq and name 
it qi on which the ground state is Qi. If qi = qo then 
qo is obviously the point we want. Otherwise we find the 
intersection point of Pi, P2 and Qi and name it q2. The 
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process can be repeated and the sequence {qn} will even- 
tually converge to the adjacent degenerate point due to 
the convexity of the ground state surface. The process of 
finding adjacent degenerate points is iterated recursively 
until it reaches the outside of the predefined region, or 
it finds a point that has already been visited. By con- 
necting degenerate points with straight lines, all ground 
states within the region are identified. 

Using the method described above, a ground state pic- 
ture on the H — A plane of a particular 32'^ realization 
(seed 1003) was computed and is shown in Fig. [3a). Co- 
existence lines are drawn with thickness reflecting the 
jump in magnetization to visualize the size of discontinu- 
ity. Most of degenerate points are inters ection points of 
four ground states, as illustrated in Fig. 6(b) The four 
ground states differ by the orientation of two separate 
domains, which are typically small, as is the discontinu- 
ity in physical quantities between them. More interesting 
are triple points where three ground states are degener- 
ate. Here a single coexistence line bifurcates into two 
lines in a Y shape, as illustrated in Fig. 6(a) The state 
on the top of the Y results from the break-up of a rel- 
atively large domain while this domain flips as a whole 
across the vertical line of the Y. A triple point has some 
characteristics of a thermal first order transition where 
two ordered state coexist with a disordered state. 

There are several thousand lines in the ground state 
picture in Fig. |5(a)| but most of these lines have small 
jumps in bond energy and magnetization. We believe 
that only relatively large jumps contribute to the singu- 
larity and, to emphasize these jumps, we simplify the pic- 
ture by removing the lines representing small jumps. In 
Fig. 5(b) is the same picture as Fig. 5(a) but a large num- 
ber of lines with small bond energy jumps {5e < 0.03) 
have been eliminated. This simplified picture reveals a 
tree-like structure built from triple points. The first or- 
der line separating the two ordered states is the trunk 
of the tree, which bifurcates at the finite size critical 
point, located at the center of the picture, into two main 
branches. Above the finite size critical point the ground 
states are disordered. The points labeled 1 and 2 cor- 
respond to the large jumps with the same labels in Fig. 
121 The inset in Fig. 5(b)| shows the details of the finite 
size critical point and two other triple ponts immediately 
above and below it. In this paper the finite size criti- 
cal point is identified as the degenerate point that maxi- 
mizes the discontinuity in the bond energy, measured by 
Se* — (|e+ — eo| + |e_ — eo|)/2, where e+, e_ are the 
bond energies of the two ordered states, and eo is the 
bond energy of the disordered state, respectively. 

We propose that the critical region of the ground state 
picture can be rescaled in such a way that pictures for 
various system sizes are statistically indistinguishable 
from one another. The required scaling involves the 
width of the pictured region Wh in the H direction, 
height Wa in the A direction and lower cut-off Scmin 
for coexistence lines retained in the picture. The scaling 
of Semin should follow Eq. (|10|l . The picture will include 
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FIG. 5: Ground states of a given realization (seed 1003) with 
system size 32"^ in the H — A plane, (a) All the ground states 
of a single realization with L — 32. The lines are coexistence 
lines of two ground states. The thickness of a line is propor- 
tional to the magnetization jump across the line.(b) The same 
realization as in (a), but only coexistence lines with bond en- 
ergy jumps Se > 0.03 are shown. Numbers 1 and 2 correspond 
to the two largest jumps shown in Fig. |21 The inset in (b) is 
a blow up of the region around the triple point identified as 
the finite size critical point and also showing the triple points 
immediately above and below it. 
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(a) 




(b) 



FIG. 6: (a) Three states are degenerate at a triple point. The 
"+" and "— " sign are used to indicate the direction of spins 
in a domain. The spins in the domain are all pointing up 

(denoted as "++") or all down (" ") in the ground states 

below the triple point, while the domain breaks up ("+—") 
in the ground state on top of the triple point, (b) Four states 
separated by two intersecting straight lines. The four states 
differ from each other in two separate domains. 



0.2 

[|dH/dA|] 




4 5 6 7 8 9 10 



20 30 



FIG. 7: Scaling of the average inverse slope of coexistence 
lines near the finite size critical point (except for the first order 
line). The slope of the best-fit line is —0.79 ± 0.04, which is 
in agreement with the predicted value (a + /3 — l)/^- 



a scale invariant part of the critical region if Wa scales 
as A — Ar, 



(15) 



The scaling of Wh is expected to be the same as the 
scaling of the external field H, which has been given by 
Bray and Moore in their scaling theory of the RFIM^, 



(16) 



In Fig.|51the parameters of the pictures are scaled such 
that Je^^L^i-")/", WaL'^^" and VF//L(2-"-/3)/'^ are aU 
held constant. Although different realizations have quite 
different ground state patterns there is no apparent way 
to distinguish between different system sizes. 

In order to test the scaling of the ground state pictures 
more quantitatively we measure [|d-ff/dA|], the average 
of the absolute value of the inverse slope of coexistence 
lines near criticality (except for the first order line) as a 
function of system size. The result is shown in Fig. |7| 
We measure the inverse slope of coexistence lines rather 
than the slope itself, because in some realizations the 
slope is very large, and thus the average of dA/dH is not 
well-behaved. The slope of the best-fit line is —0.79 ± 
0.04. From Eq. ^ and Eq. ^ we expect [\dH/dA\] - 
2^(a+/3-i)/^^ The measured value -0.79 ± 0.04 is close to 
{a + P — if (1 — q)Iv 0.84 as we have calculated 
in Sec. ^ and /3 « as generally accepted. 

We then measure the strength of the external field at 
the finite size critical point [|i?c|], which should have the 
same scaling as Wui and show the result in Fig. |S1 The 
slope of the best-fit line is —1.60 ± 0.06, which is again 
consistent with the expected value of (a + /3 — 2)/;/, if 
/3 sa 0, and our measured values of exponents (1 — aj/u w 
0.84, and l/u k, 0.8 are used. 



[|HJ] 




FIG. 8: Scaling of the average strength of the extenal field 
at the finite size critical point. The slope of the best-fit line is 
— 1.60 ± 0.06, which is in agreement with the predicted value 

(a + /3-2)//.. 



IV. POSITIVE TEMPERATURE RESULTS 

We have studied the RFIM at fixed disorder strength 
of Ao < Ac and zero external field for all T > us- 
ing the Wang-Landau algorithmS^. The Wang-Landau 
algorithm is a fiat histogram Monte Carlo method that 
also automatically determines the density of states g{E). 
Thermodynamic quantities related to energy, such as the 
specific heat, can then be derived from the density of 
states at all temperatures. In order to get the magneti- 
zation or susceptibility, we collected joint magnetization 
and energy statistics. The algorithm smooths the energy 
landscape and improves on the performance of the con- 
ventional Metropolis algorithm. Using the method we 
can determine the specific heat and susceptibility over a 
broad temperature range for systems up to size 32'^. After 
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FIG. 9: Ground state pictures for different system sizes plot- 
ted with scaled coordinates and lower bound for bond energy 
jumps. Each figure shows the scaled ground state picture for 
a single realization. 



we obtain the density of states, we use the MetropoUs al- 
gorithm to obtain average spin configurations for selected 
temperatures. 

Our first observation is that for some large enough sys- 
tems (> 16^) and strong enough disorder, the specific 
heat and the susceptibility display one or more sharp 
peaks, as illustrated in Fig. ^| For a given realization, 
the sharp peaks in the specific heat and the susceptibility 
occur at the same temperatures. The sharp-peaked tran- 
sitions have some first-order-like properties. For exam- 
ple, the energy probability density = e~^/'^ g{E)/Z 
displays double peaks, and the Binder cumulant B{T) ~ 
1 — {rti^) /^{m^) is negative at the temperature of the 
sharp peaks. The angular bracket stands for a thermal 
average. The double peaked energy distribution and neg- 
ative Binder cumulant are shown for a IB'^ system (seed 
1013) in Fig. ^] These first-order-like features result 



TABLE V: Number of realizations that have a double-peaked 
energy probability densities at their specific heat peaks (Ndp), 
and total number of realizations simulated (Ntot) as a function 
of disorder strength Ao and system size L. 



L Ao Ndp Ntot Ndp /Ntot 

8 1.5 256 0% 

8 2.0 64 0% 

16 1.5 6 96 6.25% 

16 2.0 21 96 21.8% 

32 1.5 3 9 33.3% 

32 2.0 6 9 66.6% 



from the coexistence of two states that differ by flipping 
a large domain as we will see more clearly later. Prelimi- 
nary statistics from a small sample of realizations suggest 
that the fraction of realizations showing sharp peaks in- 
creases with the system size and the strength of disorder, 
as shown in Tabled Here we call a transition "sharp" if 
the sampling probability has two peaks at the transition 
temperature. 

The sharp peaks occur at different temperatures with 
different height for different realizations, and they are 
smoothed out by an average over realizations. We show 
in Fig. ^]the average specific heat for all of the 16'^ sys- 
tems we have simulated at Aq = 2.0. Though there are 
21 sharp-peaked realizations out of a total of 96 simu- 
lated (see Table 01, the average specific heat is a smooth 
curve. The difference between the average specific heat 
and that of individual realizations shows that there is no 
self-averaging at positive temperature, similar to what 
we have already seen at zero temperature. The lack of 
self-averaging has also been observed in the bimodal dis- 
tribution RFIM^a. 



V. RELATION BETWEEN GROUND STATES 
AND THERMAL STATES 

The zero temperature fixed point picture of the RFIM 
phase transition predicts that the behavior in the critical 
region at positive temperature is determined by the com- 
petition between couplings and random fields with ther- 
mal fluctuations serving only to renormalizc the strength 
of these quantities. The results presented in this section 
suggest that a strong version of the zerio temperature 
scenario holds for individual realizations of normalized 
random fields. We will show that the sharp peaks in the 
thermodynamic quantities can be matched one to one 
with the large jumps at zero temperature. Furthermore, 
the spin configurations on either side of the sharp peaks 
can be mapped onto the ground states on either side of 
the corresponding large jumps. Similar correlations be- 
tween ground states and thermal states were found in one 
dimension^ ^. 

We illustrate the above statement for one 32'^ realiza- 
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(a) 16^, Ao = 2.0, seed; 
1013 



(b) 16^, Ao = 2.0, seed; 
1013 



^ 0.6 
§0.4 
0.2 
0.0 




-11000 -10000 
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(a) 16^, Ao = 2.0, seed; 
1013 



(b) 16^, Ao = 2.0, seed; 
1013 



(c) 16^, Ao = 2.0, seed; 
1050 




(d) 16^, Ao = 2.0, seed; 
1050 



FIG. 11: The first-order-like properties of sliarp peaks, (a) 
sliows tlie double-peaked sampl ing pro bability p{E) at tiie 
sfiarp peak for tlie system in Fig. |10(a)| (b) sfiows the Binder 
cumulant B as a function of temperature for the same system, 
which becomes negative at the sharp-peaked transition. 




(e) 32^, Ao = 1.5, seed; 
1000 



(f) 32^, Ao = 1.5, seed; 
1000 



FIG. 12: The average specific heat of 96 realizations of size 16'^ 
and disorder Ao = 2.0. Although some of these realizations 
have sharp peaks, the averaged specific heat is smooth. 




(g) 32^, Ao = 2.0, seed; 
1003 




(h) 32^, Ao = 2.0, seed; 
1003 



FIG. 10: The specific heat C and the susceptibility x of four 
realizations of the RFIM. The sharp peaks in the specific heat 
and the susceptibility of a given realization occur at the same 
temperatures. 



tion (Ao = 2.0, seed 1003) whose sp ecific h eat and sus- 
ceptibility are shown in Fig. 10(g) and 10(h) respectively. 
There are two major peaks in the specific heat and the 
susceptibility, and each of them are related to the two 
major jumps in the bond energy and the magnetization 
at zero temperature, as shown in Figs.[21and 5(b) (labeled 
as 1 and 2). 

The connection between the zero temperature transi- 
tions and positive temperature transitions is confirmed 
by the correlation between the average spin configura- 
tions near the positive temperature transition and the 
ground states near the zero temperature transition. For 
a single realization of random fields, we obtain the ther- 
mally averaged spin configuration at a given temperature 
near the peaks using the Metropolis algorithm. We start 
our simulation from the ground state, and then employ 
the Wang-Landau algorithm without modifying the al- 
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ready obtained density of states, until the microstate falls 
into an energy bin that has a significant sampling proba- 
bility at the given temperature. The Metropolis updates 
are then used to obtain the averaged spin configuration. 

Figures 13(d)| 13(e) and 13(f) show one plane through 
the system with Aq = 2.0 and at temperatures just be- 
fore peak 1 (T = 2.2), just after peak 1 (T = 2.5), and 
just after peak 2 (T = 2.8), respectively. The differ- 
ence between the states shows that the sharp peak cor- 
responds to flipping a relatively large domain. It is evi- 
dent that these three states are strongly correlated with 
the ground state spin configuration before the jump 1 
(A = 2.36), just after jump 1 (A = 2. 41), and just a fter 
jump 2 (A = 2.54), as shown in Fig. 13(a) 13(b) and 
|13(c)| respectiv ely. (T he labels of jumps and peaks are 
given in Figs. |21 |5(b)| and 10(g) ) 

Some correlation between ground states and thermal 
states persists to much smaller values of Aq in a regime 
where the thermodynamic properties no longer display 



sharp peaks. Figure [T3 (g)| [l3 (hj and 13 (i) show the same 
realization of disorder ana same plane through the 
system but with Aq = 0.5. Here the specific heat has 
a rounded peak at T = 4.375. Figures ICT g). (h) and 
(i) correspond to temperatures 4.0, 4.3 and 4.45, respec- 
tively. Although there is considerable thermal "blurring" 
in these pictures, evidence of the ground states is unmis- 
takable. 

A quantitative characterization of the correlation be- 
tween ground states and thermal states can be obtained 
from the correlation measure, 




(a) 



(b) 



(c) 




(d) 



(e) 



(f) 




(g) 



(h) 



(i) 



(A) = ;^^[sgn((s.|A,0)(..lAo,r*))] 



(17) 



where the square brackets are an average over realiza- 
tions of disorder and {si\A,T) is the thermal average of 
the spin at the ith site at disorder A and temperature 
T or, if T = 0, it is the ground state spin value. For 
each realization, the temperature T* = Tmax -1-0.1 where 
Tmax is the temperature of the maximum of the specific 
heat. Thus, for each realization, we pick a thermal state 
just above the transition temperature. Figure ITU shows 
q vs. A for sizes 16^ and 32^ and Aq = 1.5, with 96 
realizations for size 16^ and 9 for size 32^. A peak in 
the correlation occurs at A r:^ 2.65 where q « 0.75. The 
value, A « 2.65, is about 0.15 larger than the average 
A at the largest discontinuity in the bond energy for 
system size 32^. The inset in Fig. 1141 shows the aver- 
age correlation between thermal states of one realization 
and ground states of another for size 16^, which is nearly 
zero as expected. A second measure, q* is obtained by 
choosing the value A* in Eq. Ijl?!) for each ground state 
realization to give the largest correlation to the thermal 
state at T* and then averaging over realizations. We 
find that for size 32^, q* = 0.80 ± 0.06 for Aq = 1.5 
and q* = 0.85 ±0.05 for Aq = 2.0. Together, these resuh 
provide quantitative confirmation that the thermal states 
at temperatures slightly above the thermal critical point 



FIG. 13: Spin configurations near the critical points at zero 
temperature and finite temperatures for a single realization 
of normalized random fields. Each panel is the same plane 
of a 32'' realization with black representing spin down; white, 
spin up; and shades of gray, the thermally averaged spin state. 
From left to right in the top two rows, panels are at A (T) 
befor e, betw een and after jumps (peaks) 1 and 2 in Fig. |5(b)| 
(Fig. |10(g)] i. Specifically, panels (a), (b) and (c) are ground 
states at A = 2.36, 2.41 and 2.54, respectively. Panels (d), (e) 
and (f) are at A = 2.0 and T = 2.2, 2.5 and 2.8, respectively. 
Panels (g), (h) and (i) are at A = 0.5 and temperatures 4.0, 
4.3 and 4.45, near the peak in the specific heat at T = 4.375. 



are strongly correlated with the ground states at disor- 
der strength slightly higher than the zero temperature 
critical point. 

The correlation between thermal states and ground 
states near the transition is consistent with, but not pre- 
dicted by the zero temperature fixed point hypothesis^^.. 
This hypothesis predicts that the renormalization group 
fiow is to a zero temperature fixed point so that the zero 
temperature and positive temperature transitions are in 
the same universality class. However, it does not predict 
anything about the spin configurations along the critical 
line for individual realizations of disorder. If the corre- 
lation of spin configurations along the critical line that 
we observe for small systems persists to large systems. 
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A 



FIG. 14: Disorder averaged correlation g of a thermal state 
just above the transition temperature at Ao = 1.5 to ground 
states at disorder strength A for the same realization of ran- 
dom fields. Solid squares for size 16^ and open circles for size 
32^. Only a few error bars are drawn to make the figure eas- 
ier to read. The inset shows the correlation of thermal states 
with ground states of a different random field realization. 

it will support the following strong version of the zero 
temperature fixed point scenario: for a given realization 
of normaUzed random fields, the sequence of states near 
the zero temperature critical point obtained by varying 
A for T = can be mapped onto the sequence of thermal 
states near the critical point obtained by varying T for 
fixed values of Aq, Aq < Ac- 

VI. SUMMARY 

In this paper we have numerically studied the RFIM at 
zero temperature and positive temperature. At zero tem- 
perature we have extracted critical exponents from the 
finite size scaling of the several largest jumps in the bond 



energy. Our measured value of exponents (except v) are 
mostly consistent with previous values but have better 
accuracy. We have found that the heat capacity expo- 
nents a is near zero. We have also portrayed all ground 
states within a small critical region on the H — A plane 
for up to 32'^ systems. The ground state pictures shows a 
tree-like structure if small jumps are removed. Although 
the ground state pictures are not self-averaging, they sat- 
isfy statistical scaling relations. That is, within a scaled 
region in H — A plane, with scaled lower limit of bond en- 
ergy jumps chosen, the ground state pictures of difii'erent 
system sizes arc statistically similar. 

We have used the Wang-Landau algorithm to study 
the RFIM at positive temperature. This algorithm en- 
abled us to obtain the density of states and to derive 
the specific heat and susceptibility over a broad range 
of temperatures for systems up to size 32'^. We have 
observed that for some disorder realizations the transi- 
tion is characterized by sharp peaks in the specific heat 
and the susceptibility. The sharp-peaked transition has 
some first-order-like features and the fraction of realiza- 
tions that have sharp peaks increases as the system size 
or the strength of disorder increases. The sharp peaks in 
the thermodynamic functions result from flipping a large 
domain and are related to large jumps in bond energy 
and magnetization at zero temperature. More specifi- 
cally, the thermal average spin configurations near the fi- 
nite temperature transition arc correlated to the ground 
states near some corresponding large jump at zero tem- 
perature. This phenomenon suggests a strong version of 
the zero temperature fixed point scenario. It remains to 
be seen whether the correlation between critical ground 
states and thermal states persists to large systems. 
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